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O ■ ABSTRACT 

A description of the dynamics of a coUisionless, self-gravitating fluid is developed and applied to 
follow the development of Large Scale Structures in the Universe. Such description takes on one of the 
\ assumptions of the Adhesion Approximation (AA) model, i.e., the introduction of an artificial viscosity 

v in the Euler equation, but extends it to deeper non-linear stages, where the extrapolation of the 
linear relation tp = —cf) between the velocity and the gravitational potentials - at the basis of both the 
Zel'dovich and the Adhesion models - is no longer valid. 

This is achieved by expanding the relation bewteen and ^p, in general non explicitly computable, in 
powers of the small viscosity v. In this case, the evolution of the velocity potential is described by a 
■ diffusion-like equation for the "expotential" field ^ — exp{—^/2v). Such equation includes a source term 

^ ' ^(C) which expresses the relation between </) and ^ as a series expansion in powers of v. Such term is 

OO \ related to the onset of the non-linear evolution of the velocity potential and grows from zero (the limit 

04 ■ corresponding to the Adhesion Approximation) with increasing time. For terms in V{£P) up to order o{y) 

(the only ones which can be expressed in a fully Eulerian form), the diffusion equation is solved using 
the path-integral approach. 

The Adhesion Approximation is then recovered as a "free-particle" theory, (corresponding to V{$P) = 0) 
where the dynamics is determined by the initial value of ^. Our inclusion of the lowest-order term in V{^) 
substantially changes the dynamics, so that the velocity potential at a given time in a given Eulerian 
i-C ' position depends on the values taken at all previous times in all other coordinates. This is expected in 

the non-linear regime where perturbations no longer evolve independently, but "feel" the changes of the 
Q ■ surrounding density field. 

^ ' The path-integral solution is computed numerically through an algorithm based on Monte Carlo re- 

, alizations of random walks in the Eulerian space. In particular, the solution ^ at any cosmic time is 

' obtained upon averaging the value of the potential V{SP) at previous times in the Eulerian locations 

reached by the random walks. 

The solution is applied to the cosmological evolution of a Cold Dark Matter density field, and the 
' results are compared to the outcomes of an N- body simulation with the same initial condition. The 

\ velocity field in the presented extended Adhesion (EA) description is obtained numerically using the 

random walks algorithm described above. For case of a null potential V(£,) = 0, this constitutes a novel 
implementation of the Adhesion Approximation which is free from the numerical errors affecting finite- 
difference solution schemes for partial differential equation and is faster than the Gaussian convolution 
algorithm adopted by Weinberg & Gunn (1990a). When the first order term of the potential V{^) is 
included, the proposed extension of the Adhesion approach provides a better description of small-scale, 
deeply non-linear regions, as it is quantitatively shown by the computation of some statistical indicators. 
At larger scales, the satisfactory description of the large scale texture and of the voids given by the 
canonical Adhesion Approximation is preserved in the extended model. 

Subject headings: cosmology: large scale structure of the Universe - cosmology: theory - galaxies: 
formation - hydrodynamics 
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1. INTRODUCTION 

The formation of cosmic structures in the Universe is one of the key problems in cosmology. In the current, stan- 
dard theory (sec. c.g, Peebles 1993), they develop from the amplification of initially small density perturbations due to 
gravitational instability. For a coUisionless self-gravitating fluid (similar to the non-baryonic dark matter, believed to 
dynamically dominate the Universe), the description of such process is given by the usual hydrodynamical equations for 
the velocity and the density fields (Euler and continuity equations) in expanding coordinates plus the Poisson equation 
coupling the gravitational potential (j) to the density field 5. In the linear stage, when the growing density is sufficiently 
small, such equations can be solved analytically using a first order perturbation technique for the Eulcrian density and 
velocity fields; higher order perturbation theory can give insights on the quasi-linear regime where, although small, 5 
becomes comparable to unity (see, e.g., Munshi, Sahni & Starobinsky 1994); however, the non-linear regime when \5\ > 1 
is usually followed in detail through numerical N-body simulations (see Bertschinger 1998 for a review). These show that 
gravitational instability of dark matter can indeed lead to the kind of large scale structures observed in the Universe 
(starting from de Lapparcnt, GcUcr & Huchra 1986: sec Maddox ct a. 1990; Saunders ct al. 1991; Vogcley ct al. 1992; 
Schectman et al. 1996), which include filamentary overdensities and two-dimensional sheet-like structures, the so called 
pancakes; at the confluence of such one- or two-dimensional structures, high constrast (S > 200), knotty structures appear 
to form, corresponding to galaxy clusters. 

At the same time, several approximation schemes and semi-analytical approaches have been developed for studying 
the formation of Large Scale Structures (hereafter LSS). The aim of such approaches is to gain an insight into the 
physical processes leading to structure formation, to comprehend and check the results of the simulations and to provide 
a computational tool which is usually faster and easier to implement. 

The first step in an analytic approach to LSS was taken by Zel'dovich (1970), who proposed to extrapolate the linear 
behaviour of the velocity field in the non-linear regime and to express the evolution of each particle in terms of its 
Lagrangian coordinates. In this case, the trajectory of each particle evolves along the free-flight path determined by the 
initial Lagrangian velocity (for a review, see e.g., Shandarin & Zel'dovich 1989). In terms of the velocity potential (as 
long as no orbit crossing occurs the flow is irrotational and the velocity can be expressed as Vip), it can be easily shown 
that this corresponds to assume the velocity potential equal to the gravitational potential, a condition which holds in 
linear theory. The Zel'dovich approximation (ZA) turned out to work suprisingly well, even beyond the regime where its 
approximations are realistic (see, e.g., Melott et al. 1983; Shandarin & Zel'dovich 1989; Sahni & Coles 1995). However, it 
presents some shortcomings, which can be traced back to the neglection of the back reaction of the evolving density field 
on the gravitational field, and hence on the particle velocities. Thus, small-scale variations are transferred to much larger 
scales, resulting in a poor description of ovcrdcnsc; rc;gions; the collapse time of condcinsations is in general overestimated; 
the absence of a restoring force results in the formation of multi-stream regions, due to the crossing of particle orbits, so 
that pancakes indefinitely thicken after their formation. 

To overcome the above problems several improvements have been proposed. The Lagrangian perturbation theory 
developed by several authors (see Moutarde et al. 1991; Buchert & Ehlers 1993; Bouchct ct al. 1995) includes the ZA 
as a first order solution; at higher orders, the particle displacement field is determined not only by the initial Lagrangian 
velocity (like in ZA) but also by the acceleration field. When compared with N-body simulations (see Bouchet et al. 1995; 
Melott, Buchert fc Weiss 1995) the higher-orders Lagrangian approach turns out to improve the ZA for what concerns the 
first collapsing objects, the statistics of overdense regions and the compactness of clusters; on the other hand, imderdense 
regions are better described in ZA. In any case, the thickening of pancakes after shell- crossing remains as a major 
drawback of such elaborate models. 

Alternative approaches aiming at overcoming the thickening of pancakes typical of the ZA are constituted by the frozen- 
flow approximation (Matarrese et al. 1992), the linear potential approximation (Bagla & Padmanabhan 1994; Brainerd, 
Scherrer & Villumsen) and the Adhesion Approximation (AA hereafter; Gurbatov & Saichev 1984; Gurbatov, Saichev 
& Shandarin 1985; Gurbatov, Saichev & Shandarin 1989; for reviews see Shandarin & Zel'dovich 1989; Vergassola et al. 
1994; Sahni & Coles 1995). 

The flrst two propose to "freeze" the initial velocity and potential field, respectively, to their initial value; particles 
then move with a velocity determined by the local Eulerian value of the initial velocity potential (or following the line 
of force of the initial gravitational potential in the linear-potential approximation). Such approaches avoid the shell- 
crossing occurring in ZA but, as shown by comparison with aimed N-body simulations (Sathyaprakash et a. 1995), break 
down relatively early, soon after the non-linear length scale exceeds the mean distance between peaks of the gravitational 
potential; in particular, the frozen-flow approximation, although reproducing reasonably well the density probability 
distribution of the dark matter field, fails in moving the mass particles to the right places when compared with the 
N-body simulations. 

The AA, on the other hand, takes on the basic assumption of the ZA (i.e., the equality of the gravitational and 
velocity potentials) and introduces an artificial viscosity into the Euler equation to avoid orbit crossing. Although 
introduced phenomenologically, later investigations by Buchert & Dominguez (1998) show that it is indeed possible to 
obtain viscosity- like terms from kinetic theory of self-gravitating coUisionless systems (although the corresponding multi- 
stream forces are, in general, anisotropic, unlike the assumption of the AA). The effect of the viscosity ly in the limit ;y ^ 
can be straightforwardly computed; particles initially follow their linear trajectories (the same of the ZA), but when flow 
lines intersect, the colliding particles stick to each other, thus binding collapsed structures and flxing the principal failure 
of the ZA. The networks of structures resulting from implementations of the AA have a remarkable resemblence with 
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those emerging from N-body simulations; indeed, rms density fluctuations agree to better than 20 % on scales larger than 
~ 5 Mpc (Weinberg & Gunn 1990a; Kofman et al. 1992; see also Melott, Shandarin & Weinberg 1994). Such agreement 
makes AA a reliable tool for sevaral astrophysical applications concerning LSS (see, e.g., Nusser & Dekel 1990; Weinberg 
& Gunn 1990b). 

Despite of the the successes listed above in reproducing the texture of LSS, the A A model is a much less satisfactory 
description when structures at smaller scales are considered (Weinberg & Gunn 1990a; Kofman et al. 1992). The density 
field is less clumpy than appears in N-body simulations, where walls and filaments fragment into dense clumps, at variance 
with the outcomes of AA. This is of course a conseguence of neglecting the back-reaction of the particle distribution on 
the evolution of the velocity fleld. 

Here it is proposed an Eulerian approach to extend the AA to deeper, non-linear regimes. It allows both to avoid the 
shell-crossing problem of the ZA and to go beyond the approximation (valid in the linear regime) stating the equality 
between the velocity and the gravitational potential, which is the basic ansatz of both the ZA and the AA descriptions. 
As a result, the velocity field felt by a particle at any given time is affected by the dynamical evolution occurred up to 
the considered time, as it must be in the non- linear regime (see Coles & Chiang 2000). 

The equation governing the evolution of the velocity field is found by expanding the relation bewteen and ip^ in 
general non explicitly computable, in powers of the small viscosity v. In this case, the Bernoulli equation which governs 
the evolution of the velocity potential in AA can be recast (after the Hopf-Cole transformation ^ = exp{—il} /2v)) as 
a diffusion-like equation with a source term V{C) constituted by an expansion in powers of such term expresses the 
departure of the velocity field -0 from the linear behaviour -0 = —4>, assumed to hold in AA. Since the lowest order term 
in V(^) can be expressed in a completely Eulerian form, it is possible to solve such equation in the Eulerian space in this 
restricted case. 

The solution is obtained using the formalism of brownian motion, equivalent to the path-integral formulation used in 
quantum mechanics and in statistical physics. In the limit of small times, V{£^) — > 0, (corresponding to a null potential in 
the language of path-integrals) one recovers the standard AA, which thus constitutes, in the language of path-integrals, 
a free-particle theory (see also Jones 1999) whose solution is determined in terms of the initial field; the first order term 
of the potential V{C) ~ corresponding to a "theory with interactions" in the language of path-integrals - introduces the 
non-linear corrections to AA; the solution at a generic time depends not only on the initial field, but also on its values at 
all previous times and at all other coordinates. 

To test the proposed description, the solution is applied to the cosmological evolution of a Cold Dark Matter density 
field. The corresponding velocity field is obtained numerically at any time by constructing - for each Eulerian coordinate 
- a set of random trees, which are used to compute the path-integrals with specific forms of the interaction potential. 
For a null potential, this constitutes a novel implementation of the AA which is free from the numerical errors affecting 
finite-difference solution schemes for partial differential equation and it is faster than the Gaussian convolution algorithm 
adopted by Weinberg & Gunn (1990a); when the first order term in the path-integral potential is included, the evolved 
field gives a better description of small-scale, deeply non-linear regions, as shown by the comparison with an aimed N-body 
simulation. 

The plan of the paper is as follows. The Bernoulli equation for the velocity potential '0 typical of the AA model is 
introduced and extended beyond the linear evolution (Sect. 2). After a canonical change of variables (the Hopf-Cole 
transformation) such an equation is transformed into a diffusion-like equation for if) with a potential term which expresses 
the non-linear evolution of the transformed velocity potential. The latter term is expanded in powers of the artificial 
viscosity; the solution for the leading order is obtained using the path-integral formalism (Sect. 3), and it is numerically 
implemented through the construction of random walks for the transformed velocity field (Sect. 4). The results (Sect. 5) 
are then compared with the outputs of an N-Body simulation with a Cold Dark Matter power spectrum for the initial 
density perturbation field. Sects. 6 and 7 are devoted to conclusions and discussion . 

2. BASIC DYNAMICS 

It can be easily shown that the Euler equation for a coUisionless self-gravitating fluid in the Newtonian limit in the 
expanding Universe can be conveniently rewritten as a Bernoulli equation for the velocity potential, when rescaled variables 
are used (Gurbatov, Saichev, & Shandarin 1989; Kofman 1991). If the comoving peculiar velocity field u = Viia is 
expressed in terms of the gradient of a velocity potential ip and of the time derivative of the expansion factor a, then the 
evolution of the velocity potential is governed by 

|^ + i(vv,f.-JL(. + „, (1) 

Here ^ is the gravitational potential divided by 3to/2aQ, where to and ao are the initial time and expansion factors, 
respectively. [| 

The ZA can be recovered from eq. (^ imposing that (p = ""0? ^ condition which is valid in the linear regime. The 

^ Whenever irrelevant for the exposition and for the computations, the spatial argument x will be omitted; the dependence on the expansion 
factor will be often indicated with a subscript, so that i/)(x, a) will be often indicated as ^a- The subscript will be used for fields computed 
at the initial time, so that ijiQ = ipag ■ 
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solution of eq. (|^) is then 



V'(x,a) = ?Ao(q) 



(x- 



(2) 



2 (a - ao) 

where x and q are the Eulerian and Lagrangian coordinates of a particle with trajectory x(q, a). The solution of eq. (0) 
is characterized by the remarkable property V^iplx, a) = VqV'o(q)- Thus the particles trajectories are flee-flights with 
X = q + (a — ao)ug determined by the initial velocity = VV'q = —V(f>q. 

The AA approximation consists in keeping the ansatz ip = —ip, but adding to the l.h.s. of eq. (|l|) a viscosity term 
— With the Hopf-Cole transformation ip = —2vln£^, the Bernoulli eq. (0) is transformed into a linear diffusion 
equation d£,/da = vV^^ where the expansion factor plays the role of time; the solution is well known to be the convolution 
of the initial condition with a Gaussian whose variance is proportional to the time variable. Transforming back the solution 
for ^ into the velocity potential one obtains the expression for in AA, which reads 



ip{-x, a) 



-2i> In 



(Tqe 



5(x,q,a) 



where the action is 



S'(x,q,a) = %po{q) + 



(x-q)2 



2(a - ao) 

It can be shown (see, e.g., Vergassola et al. 1994) that, in the limit ^ 0, the solution reads 



(3) 
(4) 



a'0(x, a) = sup 



aiJai^) -q^/2 + x-q 



x'12 



(5) 



The confluence of different Lagrangian points into a single Eulerian coordinate gives rise to the formation of caustics 
and knots, reproducing the skeleton of LSS, and avoiding the shell-crossing. Indeed, given a coordinate x at a time 
corresponding to a, the Lagrangian points corresponding to orbits leading to x are all the q* where the maximum in eq. 
(||) is attained. Shell-crossing does not occur in AA since the property [q*(x, i) — q*(x',t)] • (x — x') > holds. 

3. EXTENDING THE ADHESION APPROXIMATION 

All the schemes discussed above are characterized by the extrapolation of the relation — -~(p to the non-linear regimes, 
with the resulting limitations discussed in the Introduction. 

A step forward can be made considering the remaining equations for the dark matter fluid. The continuity equation 
can be recast in terms of the rescaled density field 77 = 5 -I- 1 to read 



dri 

— ^ u • V77 + 77V • u = 0. 

oa 



(6) 



A formal solution of the above equation can be found upon integrating along the particle trajectory x(a), where x is the 
comoving Eulerian coordinate. Then one obtains 



5(x,a) = [So{q) + l]e 



■ I da' V-u{x(q,a'),a') 



1 



(7) 



where integration over Ca{x) indicates integration over the particle trajectory from the Lagrangian coordinate q at the 
initial time ao to the Eulerian position x at the time corresponding to a. The above density field is related to the 
gravitational potential by the Poisson equation 



(8) 



To obtain an equation for ip we start from eq. (Q) modified with the addition of the viscosity term —vS/'^ip on the l.h.s. 
After substituting to the value 5{x, a) /a obtained from eq. (p[) one finally obtains 



{do{q) + l)e •^<^»<-) 



(9) 



To transform the l.h.s. of eq. 
to obtain 



into a diffusion term, we perform the canonical Hopf-Cole transformation ip = — 2i/Zn^, 



A 



21^ 2v' 



AC 



3 



[I + 2ai^AZ<o(q)] e 



2v f 
Jc 



A In^ da 



2i/A 



(10) 



where we have used the property (5o(q) = — aA^o(q) = 2i^a Aln^o(q) valid in the linear regime, as appropriate since the 
above quantities are computed at the initial time corresponding to Cq. 
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We seek for a perturbation expansion of the r.h.s. of the above equation, which at lowest order must be zero (according 
to hnear theory and to the ansatz in ZA and AA) and at higher orders detach from the null value according to the growth 
of structures in the non-linear regime. To this aim, we expand the exponential on the r.h.s. in powers of the small viscosity 
u. Keeping terms up to 0{v'^) (the order of the diffusion term on the l.h.s.) one obtains 

^^_^A^= A^y(^) (Ua) 



V{0 ^ A- 



da 2 a 

aA Zn^o(q) + /^^^^^ A In^ da' + {v/2) ( J^^^^-, A In^ da'f 



+ 



2vt\ Inio (q) / Aln^da' - AlnS, 

JCa{x) 



(lib) 



The computation of all terms involving the inverse-Laplacian operator A ^ is of course extremely difficult. It is possible 
however to put in evidence some of the terms in the "potential" V{^). In particular, it turns out (see Appendix) that 

A-i Qa = A~^ \ [ Aln^ da'] = / da' In^a' + 2uA-^ / VQ ■ V In^a' da' (12) 

'-JCa{x) -' Jao J 

and 

A/<o(q) = A/< + 2r/ / V {A In^) ■ V In^ da' . (13) 

JCa{x) 

Inserting the above crjuations into cqs. (11), these can be written in compact form as a diffusion equation with a source 
term (or Shroedinger-like equation in Eucledian time) 

|=.Ae + A^(y,+,^,) (14a) 
Vi= f ln^a'da'/a = Jn^ (14b) 

J an 



Vo = A-1 



2JVQ-Vln^^,da' (/c„(:.) A/n^rfaQ^ 2 A lne»(q) /^^^^^ Aln^da' 
a 2a a 



+ 2 / V{AlnO -Vln^da' . (14c) 

The above equation shows that the non-linear dynamics with viscosity is modified with respect to the AA by the effect 
of two terms, i) A sort of "time average" of the velocity field, corresponding to "potential" Vi; this constitutes the 
lowest-order modification to the linear ansatz ^l) = —(p of AA. ii) An explicitly non-local term V2, involving the inverse 
Laplacian; of course, it is overwhelmingly difficult to compute any of the terms contained in V2. 

In the following, the solution of eq. (14) will be restricted to the first term Vi, the one which is treatable in a fully 
Eulerian form; we shall refer to the corresponding dynamics as Extended Adhesion (EA) model. Note that the other term 
V2 enters the equations multiplied by the small viscosity v, so that its effects on the dynamics should be small. Although 
this is an encouraging property, a rigorous full perturbative solution of eq. (14) would require keeping the terms up to 
order v since this is the order of the diffusive term vAS^ characteristic of adhesion models. 

The lowest order (in v) term of the potential (the "time-average" in Vi) is expected to introduce relevent modifications 
to the dynamics with respect to AA, thus representing a significant step forward in the description of the non-linear 
regime. First, we note that the source term V\ is growing from zero with increasing time (Vi — > for a — > Oo) as is 
expected for a term connected with the departure from linearity. Thus, for small times the system behaves much like the 
AA; at later times the term Vi will set in, affecting the dynamics in the deeper non- linear regime. The first modification 
introduced by V\ is that the velocity field at a given time is no longer determined by the initial field ^o, as was resulting 
from the linear ansatz if) = — ^ characteristic of the ZA and AA (see eqs. 2-4). 

At early times, in the quasi-linear regime, further insight into the potential term Vi can be gained by a series expansion of 
the velocity potential ijj = ijj^'^'^ +1/"^^^ around the initial time oq; terms of increasing order correspond to consider increasing 
powers of a. The dynamics resulting from eq. (14a) and (14b) can be compared to the exact perturbation theory (valid 
for small density contrasts and close to the initial time) up to 2nd order in ■0 (where exact solutions are available). The 
first order term (the linear solution) is the same for the exact theory, AA and EA, namely Tp'''^\x) = —(j)Q{x)\ the second 
order term in AA and EA can be obtained by inserting the Ist-order solution in the (VV')^ term in the Bernoulli eq. 1 (at 
early times particles sticking occurring in A A and EA is unimportant, so that the viscosity term can be neg_lected in this 
restricted context) and letting = —(j) (for AA, according to the Zel'dovich approximation) an <j) + if) = for the EA 
(after transforming eqs. 14a and 14b back into the tf) variable), where V' is the average over a of the velocity potential which 
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we express as V' ~ After the above procedure, the AA and EA solutions at 2nd order are f/'^^ = — (l/2)(V0o)^ « and 
tl^^^\ = — [/3/(3 — 2/3)](V(/)o)^ a; these have to be compared with the exact 2nd order solution which contains a non-local 
term (i.e. 3/7A~^ [V(Vc6oA(/)o)], which is not reproduced by AA and by EA with the potential Vi) plus a local term 
V'exLt = -(6/21)(V(/)o)^ a (see Munshi & Starobinsky 1994). Note that (for /3 « 1/2 at 2nd order) the EA at second order 
yields i^EA ~ ~(V4)(V(/)o)^ a which is closer to the exact term than AA. Thus, when a time expansion of the velocity 
potential is considered for small times, the potential obtained from both EA and AA have the same spatial structure than 
the local part of the exact 2nd order correction, but the EA is closer in normalization to the exact solution. 

At later times, a property of the dynamics described by eqs. (14) is that the changes of the velocity field in the course 
of the evolution now explicitly affect the dynamics. In addition, as will be discussed in detail in §4.3, the presence of a 
potential V{£,) depending on ^ in eq. (14) introduces non-local features in the solution, so that the field ^(x, o) depends 
on the value of the field at other Eulerian points. This feature is expected to arise in the non-linear regime, since the 
density fluctuations cease to evolve independently and "feel" the effect of the whole mass distribution (we refer to §7 for 
a more extended discussion on the effects of the term Vi on the overall evolution of LSS). 

To quantitatively explore the above effects we now proceed to solve the diffusion equation (14), restricting to consider 
only the first term Vi. 

4. SOLVING THE DIFFUSION EQUATION; THE RANDOM WALK APPROACH 

To discuss the solution of the diffusion equation (14) let us start with the simple case when no source term is introduced 
{Vi = V2 = in eq. (14)); this corresponds to the AA. The approach used for this case will be then extended to include 
the term Vi(^). 

4.1. Free- Diffusion: Recovering the Adhesion Approximation 

It is well known that the linear diffusion equation d^/da — uA^ = describes the time evolution of the probability 
distribution for a Gaussian random walk. Let us define a random variable bg which, as time is incremented by a step ds, 
increments its value by a random amount Sh extracted from a Gaussian distribution with variance a"^ = 2^ ds along a 
path whose time coordinate s ranges from Oq to the time a. Then the solution of eq. (14) in the point x can be written in 
terms of the initial field computed at the locations b(a), i.e., the coordinate reached by the random path by the time 
a. In particular the solution writes (see, e.g., Garter & Molchanov 1991) 

^(x,a) = (^o(b(a))) , (15) 

where the average refers to the ensemble of paths bg departing from x at time 0. We indicate with b(a) the value of the 
random walk at time a. 

Figure la illustrates how this solution works, for the one- dimensional case where a simple visualization is possible. To 

find the function ^(x, a), a number of realizations of the random walk are started from the point x at Oq. At the time a, 
each Eulerian point reached by the J — th realization bj(a) of the random walk is "projected" at the initial time Oo (the 
points bi, b„ in fig. la) and there the initial ftmction ^o(bj) is computed. The solution is obtained after averaging 
over all the possible bj, weightening with their probability to occur. Since for a Gaussian random walk this is a Gaussian 
with variance ^^(a) = 2v{a — Uo), the projections on the Uo axis of the points bj(a) (i.e., the points bi, b„ in fig. 
la) will deviate from the initial position x of the random walk with a probability 

^(b) = , ^ (x-b)V2(a-a„)_ (Ig) 

(47ri^(a - ao)3/2) 
According to what said above, the solution at time a is then 

e(x,a)= [ d^bP(h)Uh) . (17) 



Indeed, performing the Hopf-Cole transformation back to the velocity potential tjj = —^vln^, the velocity potential of 
the AA (eqs. 3, 4) is obtained. 

Note that the above solution can be written in the language of path- integrals, widely used in quantum mechanics 
(Feynman & Hibbs 1965). 

^(x, a) = J K{x, a, Xo, Uo) ^(x, a^) d{xo) . (18) 
The kernel K is the particle propagator which is generally written as 

K{x, a, x„, ao) = C e^I^Wl V[h{s)] . (19) 



The integal on the r.h.s. is actually a sum over all the random paths b{s) that connect xq at the initial time to x at the 
present time, the variable s corresponding to the time variable of the random walk. The symbol ©[^(s)] implies integration 
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over positions at intermediate times in the random walk; for a discrete walk constituted by n-steps labelled si,..., s„, 
it takes the form P[6(s,i)] ~ II^^q bsJ\/4:TT (s^+i — Si). The action 5'[6(s)] actually weights the paths contributing to 
the integral. For the free diffusion equation the action is that of a free particle, containing only the " kinetic" term 
S = — J [l/Av) [dh{s) / ds^ ds . Thus, in this language, the AA (leading to a free-diffusion equation for the transformed 
velocity field ^) corresponds to a free-particle theory. 

4.2. Diffusion with a Source Term 

The solution for a diffusion equation with a source term can be obtained generalizing the action in eq. ( p^ ) to include 
the presence of a potential term. When F(^) is included the action takes the form 





' dh{s) 


[St. 


ds 



Vi^h, s))}ds. (20) 



When inserted into the path-integral (eqs. 18-19), this provides a solution to eq. (14). Thus, the terms Vi in eq. (14b) 
constitutes an interaction potential proportional to the time average of the ijj. The non-linear effects in the evolution of 
the velocity field of a self-gravitating fluid with artificial viscosity is then mapped into a theory with interaction for the 
field ^. 

The random- walk representation of the solution defined by eq. (|2^) , is the analogous of eq. ( p^ ) and reads (see Garter 
& Molchanov f99f ) 

C(x,a) = {^o{ha)e-'-o ' ) . (2f) 

Of course, since the function ^ itself appears as an argument of the potential on the r.h.s., eq. ( pi) ) actually represents an 
equation for ^, which is equivalent to eq. (14). To show such equivalence and to discuss how theaoove solution works, let 
us write the time evolution of the field $(x,a) satisfying eq. (pi]): 

s,^,^ , - ,sov"V" . - ^^^^''^'■^^''') . (22) 

Expanding both the first and the second factor in the average at the r.h.s., one obtains 

eo(b(a + da)) = UHa)) + ■ VS,o{b{a)) + ^('^b)' Co(b(a)) (23a) 

nms))ds ^ vmis))<is^^ ^ ^^^^^^^^^ ^^sb) 



We insert the above expansions into eq. (^2|) and perform the average over the distribution function ^'(^b). If this is 
symmetric and with variance 2vda (we choose it to be a Gaussian), then the terms proportional to 5h cancel out, and 
we are left with 

^(x, a + da) = iUHa)) ^(«(''(^»'^^) + (e„(b(a)) Jl ^(«(''(^»'^^) y(C(a))rfa+ 

+ i.daV2(eo(b(a))eC^(«^(^»'^^) (24) 

out to order o(da}). After substituting eq. (Elf) for the ensemble averages, dividing by da and taking the limit c?a — > 
yields eq. (14) for a generic potential y(^) on the right-hand side. This shows that eq. ( plj ) is a reformulation of eq. (14) 
in terms of random walks. 

4.3. Implementing the Solution of the Diffusion Equation with Source Term 

Here we shall take advantage of the formulation (^ ) to develop a numerical method for computing the solution of eq. 
(14). This will allow to avoid the use of finite-difference schemes for integro-differential equations which are characterized 
by delicate numerical instabilities. 

To solve eq. (|l|) with numerical realizations of random walks we first set up a grid of three-dimensional coordinates x 
and of time steps, where the transformed velocity potential ^ has to be computed. Then we proceed through the following 
steps. 

( i ) At the initial time Oo, for each Eulerian position x, we assign the initial velocity field, and hence the initial field 
^o(x) = ^(x, 0). We initialize Nreai realizations of random walks, associated with the considered Eulerian coordinate 
X with the initial condition 6 j(ao) = x, where J is the label of the realization, J = 1, Nreai- The initial value of 
the potential Vi(^o, Oo) in eq. (14) is set equal to zero. 



( ii ) We increment the time step by da. For each coordinate x we update the random walk 6j(a) associated to it by 
extracting the increments Shj from a Gaussian distribution (with variance 2 v da) for each reahzation J. For each 
coordinate x, we update the variable &j(a) = 6j(a — da) + Sjh for each realization J of the random walk. 

( iii ) We compute £,o{b^{a)) by interpolating the initial field hi the point foj(a). 

We evaluate the action Sj{a) = Sj{a — da) + Vi(^^a-da{b^{a ~ da))^ da, entering the solution (21) using the value 
of the field ^ at the previous time step to compute the potential Vi(^). 

J) by summing up all the realizations of the random 



( iv ) We compute numerically the average definig the solution eq 
walk: 

N. 



e.(x) = -— J2 eo(fej(a))e^5(") 



(25) 



J=i 



( V ) Having found the solution at the time corresponding to a, we iterate from step (ii), until the final time is reached. 
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Fig. 1 - Panel a): An illustration of the random- walk solution to the diffusion equation with no source term; for the sake of simplicity, the 
illustration is restricted to a one-dimensional space (indicated as X in the label on the orizontal azis). To obtain the solution at the point x 
at time a, different realization of a Gaussian random walk are constructed, with the condition that they all start from the Eulerian point x at 
time Oq; in the figure, three such paths are schematically illustrated and labelled bi{a), 62(1)1 b-jio-)- The points bj(a) reached by the random 
walk at time a are then projected backwards to the initial time ao (and labelled 61, 62, f'3 in the picture). The average of the corresponding 
values of the initial field go at such points (indicated as £,o{bi), ...,^0(^3)) yields the solution. 

-Panel b): The corresponding graphical representation of the solution when a source term (of the kind of V\ in eq. (14)) is introduced in 
the diffusion equation. Here, for the sake of simplicity, only one realization (indicated by fej(a)) of the randow walk is showed. To obtain the 
solution at time a, besides computing the initial field at the location 6(a) as in Panel a), all the field £,^1 at times a' < a are required, and have 
to be computed at the points reached by the random walk at the times a' . Such values of the fields are used to compute the potential Vi, see 
point iii) in the text. Again the solution is found upon averaging over all the realizations of the random walk. 
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Analogously to the free-diffusion case, we give a graphical representation of the solution corresponding to the above 
alorithm (fig. lb). Again a random walk starting from x at a = is drawn, and the initial field is computed at 
the projected points corresponding to the Eulerian position reached by the random walks ba. However, in this case, for 
computing the field at time a the function ^a(b(a' < a)) has to be computed at all previous times (see point iii), since 
they constitute the argument of the potential Vi entering eq. (pi]). 

A very important point emerging from the above solution of eq. (14) is that the value of the field ^(x, a) at a given 
coordinate x does indeed depend on the value of the field at other coordinates. In fact, for computing the field at the 
point X, the functions ^a'<a entering the potential Vi, must be computed at all points b(a'). This is at variance with the 
case of simple diffusion, where all is needed to compute the solution is the initial field and the final localization of the 
random walk b(a). Of course, such property enters only when the "interaction" potential Vi is set in. This shows that 
even the first order term Vi in eq. (14) introduces the typical effects of the non-linear dynamics, i.e., the influence on 
the velocity field of a) the changes of it occurred at all previous times, and of b) the value taken by the field in all other 
Eulerian points. 

Thus, we expect that the resulting dynamics "feels" , at some level, the changes of the fields in the course of evolution, 
a feature which is completely missing in the AA, as discussed above. 

5. RESULTS AND COMPARISON WITH N-BODY SIMULATIONS 

To test the above method for solving eq. (14) and whether restricting to the term Vi in eq. (14) gives an accurate 
description of the non-linear dynamics, we compare the outcomes of the proposed description with those of an N-body 
cosmological simulation. Although a complete, systematic comparison between the N-body and the semi-analytic descrip- 
tions - including, e.g., different cosmological/cosmogonical initial conditions - is out of the scope of this paper, we will 
compare the density distribution and some statistical indicators which are believed to describe to some extent the matter 
field in the non-linear regime, for a given cosmological initial condition. 

To perform the simulation, we adopt an adaptive P3M N-body code (see Hockney & Eastwood 1981; Couchman 1991 
for a detailed description) for a self-gravitating, coUisionless dark matter. In particular, we use the public version of the 
Couchman's adaptive P3M code (Couchman, Thomas, & Pearce, 1995) for the evolution of the dark matter, which was 
also used to generate the initial velocity field which is evolved both by the N-body and (according to our description) 
after eq. (14). 

To emulate the behaviour of the cosmological dark matter fluid, a distribution of 64'^ particles is evolved in a comoving 
simulation box with periodic boundary conditions; the initial positions are assigned to be at the center of the cells of 
a 64"^ cubic grid. The initial displacement (velocity) is given by Uq = V'i/'o o, da (we use a definition of the velocity 
potential rescaled to the Hubble expansion, see §2). The initial velocity potential is derived under the approximation 
(valid in the early, linear regime) 4'o = —4'o, where the initial potential is a Gaussian random field with power spectrum 
■p0(fc) = A k'^~^ T'^ (k) , as it is commonly taken for primordial cosmological perturbations. The transfer function T{k) 
depends on the nature of the dark matter field; here we adopt the form appropriate for Cold Dark Matter (Davis et al. 
1985; for recent fitting forms see Eisenstein & Wu 1999 and reference therein). The spectrum is normalized to the data 
from COBE (Bunn & White 1995; Stompor, Gorsky, & Banday 1995). We assume a flat cosmology with matter density 
parameter = 1 and Hubble constant h = 0.5 (in units of 100 km/s Mpc). The physical length of the simulation cube is 
L = 6Ah~^ Mpc; the gravitational force was softened at small distance; the adopted softening parameter corresponds to 
0.2 in mesh units. 

To evolve an initial spatial distribution of particles according to our description, eq. (1 4) is solved for both the case 
of free-diffusion (corresponding to AA), and for a source term given by Vi defined by eq. (jwE) [the Extended Adhesion 
(EA) model]. At each cosmic time, the solution of eq. (14) for the transformed velocity field ^(x) is found by generating, 
for each x, a number Nreai of random walks b(a) through a Monte Carlo procedure, as described in detail in §4.3. After 
transforming back to the velocity potential ijj — —ivlnS^, the position x(a) of each particle is then updated at each time 
step to the new position x(a -|- da) = x(a) -I- Ua(x) da = x(a) + Vtpa{^) d da. 

The space grid used for the Monte Carlo solution described in §4.3 is taken to coincide with the 64'^ simulation box. 
As for the number of realizations of the random walks, it has been checked that convergence in the solutions is obtained 
already for Nreai > 10^ (the latter value requiring « 100 Mbyte of computer memory; of course, the larger is the value 
of Nreai the larger is the requested memory and the slower is the numerical implementation); performing a test for the 
free-diffusion case with Gaussian initial conditions a value Nreai = 10^ yields errors 5£,/S,exact < 10^"^ when the numerical 
solution is compared to the exact one S,exact, whose analytical form is known in this case; the results shown below are 
obtained for Nreai = 100. To numerically implement our description, we adopt a finite value of the viscosity. Since we 
adopt the adimensional expansion factor a as the "time" variable in the equations for the velocity field (see eq. 1; eq 9 
and following) the viscosity has the dimension of Length^; we adopt the value u = 10~^ pixel^ (the pixel corresponding 
to the mesh size), which ensures convergence in the sense that results with smaller values of i/ are indistinguishable; a 
discussion on the physical effects of adopting different values of v in AA is given in Weinberg & Gunn (1990a). 

While in the case of A A and EA the velocity potential i/'a is obtained from eq. (14), we recall that for N-body 
simulations the same displacement is obtained from the particle density field after solving the Poisson equation on the 
grid and integrating the resulting acceleration field. Thus, in such simulations, once the particles are moved, the resulting 
density field has to be re-computed to allow for solving the Poisson equation at the next time step. Such procedure, as well 
as double Fourier-transforms required to compute the solution of the Poisson equation, is not needed in the semi-analytic 
approaches, like AA ora EA, making them usually much faster than the simulations. In our case, the main source of 
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time-consumption in the numerical implementation is due to the large number Nreai of Monte Carlo realizations of the 
random walk needed to obtain reliable averages in eq. (pi]). 

The simulations are started at an initial time corresponding to an expansion factor Oo = 1/16 (normalized as to yield 
a = 1 at the present time). The resulting particle distribution at the final time a 1 is shown in fig. 2 for the N-body 
simulation (top panel), the AA (middle panel) and EA (bottom panel) for a slice 4 Mpc thick. The set of parameters 
adopted for the N-body simulation and for the AA end EA implementation are recalled and summarized in the caption. 
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Fig. 2 - Slices through the particle distribution at redshift z = from the N-body simulation (top panel), the Adhesion Approximation 
(middle panel) and the Extended Adhesion model (bottom panel) proposed in this paper. Each slice is 4 pixels thick, while the simulation box 
occupieas a 64^ grid. Initial conditions from a COBE-normalized CDM power spectrum in an £7 = 1 Universe with Hubble constant h = 0.5. 
The parameters used for the N-body and for the Monte Carlo implementation of AA and EA are the following: physical lenght of the simulation 
cube L = 64 Mpc; initial expansion factor ao = 1/16 (assuming the a = 1 at present time); the Plummcr softenig parameter adopted in 
the N-body simulation is 0.2; the number of Monte Carlo realizations used to obtain the velocity field for AA and EA is N^^^l = 10^; the 
artificial viscosity is f = 10~^ (mesh units)^. 

Compared to the simulation, AA reproduces well the general texture of LSS, but the small scale features are undere- 
produced. In particular, while in the simulations extended structures appear fragmented into dense knots, in AA they 
appear more as a continuous filament. This is because the effects of the changes of the field in the course of evolution 
(typical of the non-linear regime) are neglected in AA; since small scales are those which are evolved more deeply into the 
non-linear regime, its is natural that AA does not reproduces them in detail. On the other hand, EA seems to provide a 
more satisfactory description of the density field down to small scales; knotty, small-scale features are remarkably similar 
to those arising from the simulation, as is apparent, e.g., from the structures just above and below the large void at the 



11 



center of the picture. Most of the structures appearing in the simulations are reproduced by EA which seems to reproduce 
quite well the various degrees of dumpiness. 

Thus, EA seems to improve the Adhesion approach in that it provides a better description the fragmentation of filaments 
in correspondence of the denser knots. Such interpretation is confirmed by the more quantitative analysis performed in 
fig. 3, where it is shown the deviation of the density field computed in AA (top panel) or in EA (bottom panel) from that 
resulting from the N-body simulation. 




Fig. 3 - The overdensity (AT — Arjv6ody)/^iV6ocij/ of the particle density N resulting from AA (left panel) and from EA (right panel) with respect 
to the particle density from the N-body simulation Nj^i^^y is plotted as a function of the x-y coordinates, for the same z-plane of of fig. 2. 
The x-y region of the maps is a blow-up of the central region of the slice in fig. 2, smoothed with a radius of 2 pixels. 

The comparison is performed on the same slice shown in fig. 2, but limited to the region surrounding the big central 
void to provide a more clear graphical rendition. Note that the deviations of the AA density field are considerably larger 
than those occurring in the EA. Most important, while the map of the deviation in the EA shows no obvious spatial 
structure, the deviation map of AA clearly shows larger deviations correlated with the location of the filaments; even the 
perimeter of the large central void of fig. 2 can be recognized in the AA map (top panel) of fig. 3. Again this is related 
to the the lack of fragmentation of filamentary structures typical of AA. 

The above differences, of course, can be traced back to the modification of the AA velocity field induced by the 
"potential" term Vi in eq. (14). To show this in detail, the velocity field in AA (top panel) and EA (bottom panel) is 
represented in fig. 4. For better readibility, the plot refers to a further blow-up of the slice in fig. 2, namely, the region 
just above the big void (coordinates are specified in the caption), where the A A and EA yield clearly different degree 
of dumpiness. Inspection of fig. 4 shows that the main, large-scale streams arc indeed very similar; however, in EA a 
modulation of such large-scale flows appears, resulting into a break-up of the coherent motions (defining the filamentary 
regions) into more structured velocity configurations, which are responsible for the formation of knots along the filaments. 

The above consideration about the mass distribution in the different schemes can be tested more quantitatively though 
the computation of some basic statistical indicators. In particular, the correlation function S(r) and the rms density 
(5^)^/^ arc computed as a fmiction of the scale r, where the latter is obtained by counting the density of particles in 
cells with radius r, and averaging over the simulation volume. The results are shown in fig. 5 for redshift z = 1 and 
= 0. Note that, while at large scales the A A gives a fair description of the density field at r ^ 5 Mpc, at small scales it 
underestimates both the correlation fiuiction and the average density, as already obtained by Wcimberg & Gimn (1990). 
The same statistical properties seem to be well reproduced by EA, as it is shown by the agreement with the N-body 
results in fig. 5 which is preserved down to the resolution limit of the simulations. Again, this is due to the fact that 
structures in AA arise directly from features in the initial conditions, while EA, to some extent, captures the effect of the 
changes that the matter field undergoes in the course of evolution. 
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Fig. 5 - The two-point correlation function is shown in the top panels for the N-body simulations (solid line), the EA (dashed line) and A A 

(dotted line). In the bottom panels it is shown the rms fuctuations of the density smoothed with a Gaussian of radius r. The latter is expressed 
in pixels units; different lines correspond to the models as above. Left column refers to z = 1, while right column to z = 0. 
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6. CONCLUSIONS 

A description of the dynamics of a coUisionless, self-gravitating fluid has been developed and applied to follow the 
development of Large Scale Structures in the Universe. Such description takes on one of the assumptions of the Adhesion 
Approximation (AA) model, i.e., the introduction of an artificial viscosity in the Euler equation, but extends it beyond 
the approximation which make it strictly valid only in the linear regime, namely, the assumption of equality between 
the velocity and the gravitational potential, tp — —4>. The key points characterizing the proposed approach (Extended 
Adhesion, EA) can be summarized as follows. 

1) The dynamics emerging from such a novel description is determined by a diffusion- like equation for the transformed 
velocity potential ^ = exp{—'4! I2v). Such equation includes a source term ViX) (or an "interaction" term in the action, 
if the diffusion equation is considered like a Shrocdinger equation in Euclidean time) which grows from zero (the limit 
corresponding to the AA) with increasing time; this in fact describes the onset of non-linear evolution of the velocity 
potential. The AA is ten recovered, in the path-integral language, as a free-particle theory. 

2) When the "potential" V{^) is expanded in powers of the small artificial viscosity v, the term corresponding to the 
lowest order can be expressed in a fully Eulerian form. In this case it is possible to compute a solution for ^ based on the 
realization of random walks in the Eulerian space. The solution at the time a and at the point x is related to a proper 
sum over the fields computed at the preceding times at the Eulerian coordinates reached by a Gaussian random walk 
starting from x at the initial time. 

3) Such solution of the diffusion equation explicitly shows that the source term introduced in the proposed extension 
of AA affects the dynamics in two, key respects: i) the velocity potential at a given time is no longer determined by its 
initial form, but depends on the values taken at all previous times; ii) the value taken by "0 at a given point depends on 
the value taken at the other points. Both these features are characteristic of the non-linear dynamical regime, when the 
density and fluctuations cease to evolve independently and "feel" the effect of the whole mass distribution. 

Thus, the proposed extension of the adhesion approximation is expected to provide a better description of the regions 
which underwent a deeper non-linear evolution. To test such expectation, the solution for the velocity field in the extended 
Adhesion approach has been used to compute the time evolution of a cosmological dark matter field, and the results have 
been tested against N-body simulations. 

When one restricts to consider a null source term V{^) = 0, the Adhesion Approximation with finite visosity is recovered. 
In this case our results and comparison with the N-body simulations yield results similar to those in Weinberg & Gunn 
(1990a), as is shown qualitatively by the particle spatial distribution (fig. 2) and quantitatively by the correlation function 
and mean overdensity plots. However, the random walk approach adopted here yields shorter computational time than 
the Gaussian convolution method adopted by the above authors, resulting into a gain over the simulation a factor 2-3 
larger. 

When the "interaction" term Vi is set in, the evolution of the velocty field in the extended Adhesion approach is 
succesful in reproducing most of the features emerging from the N-body simulations, including the fragmentation of large 
scale structures into dense lumps. The correlation function an the mean overdensity as a function of scale resulting 
from our Extended Adhesion model agree remarkably with those from the N-body even at small scales (see fig. 5). At 
larger scales and for underdsense regions, our extension of the Zel'dovich ansatz -0 = —tp leaves invariant the satisfying 
agreement of the Zel'dovich and Adhesion approximation with the outcomes of the simulations; this is at variance with 
other semi-analytic aproaches to LSS (like the Lagrangian perturbation theory at second order, Bouchet et al. 1995). 

The above results indicate that the lowest-order term Vi in the equation for the velocity potential is effective in 
capturing some relavant features of the non-linear evolution of the velocity field. The physica l me aning of such term, as 



compared with the higher-order one V2 entering eq. (14), is straightforward. Inspection of eq. (14b) and of the Hopf-Cole 
trasformation ^ = exp{—'4) /2v) shows that such term corresponds to considering the effect of the time average of the 
velocity potential in the course of evolution. Thus, the solutions presented here correspond in a sense to "mean field" 
solutions. It is not surprising, then, that these constitute the lowest-order correction to the "free-particle" behaviour, 
corresponding to AA. The consideration of higher-order terms in V{£) would then correspond to consider the detailed 
effect of each single particle on the evolution of the Eulerian velocity potential. This effect, of course, is related to the 
detailed history of each particle and correspondingly it is expressed by non-local terms which involve integrals over the 
particle trajectory, like those in the term V2 entering eq. (14). 

7. DISCUSSION 

As recalled above, the dynamics described by the EA through the time-evolving velocity potential derived from eqs. 
(14a) and (14b) allows to improve the description of the high-density regions, more deeply evolved in the non-linear 
regime. This improvement allows i) to extend the insight on the physics of LSS to include the evolution of higher-density 
regions, and ii) to extend the cosmological applications of ZA and EA to a larger density range including overdensities 
up to (at least) 5 ~ 10 (see fig. 5) where EA (at variance with ZA and AA) still provides a satisfactory description. We 
shall now discuss the above two points in turn: 

Previous works based on N-body simulations (Melott, Weinberg & Gott 1988) and on implementations of AA (Weinberg 
& Gunn 1990a) suggested that on sufficiently large scales the process of non-linear gravitational evolution may be viewed 
as a smoothing proces on the initial density field; indeed, the results for the AA (Weinberg & Gunn 1990a) showed that 
the density field resulting at a given cosmic time is well approximated by the initial density field smoothed over a scale 
corresponding to w Aa;/3 where Ax is the average particle displacement at that time (i.e., the average comoving distance 
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that a particle has moved from its initial position). The AA allows to pin down the origin of such a behaviour; indeed, 
for large scales where AA is a satisfactory approximation, the diffusion term in the Burgers equation for the velocity field 
(or equivalently in the equation 1 for the velocity potential) has the effect that the velocity field at a point incorporates 
the contributions from the surrounding patch of initial conditions. This is explicitly shown by the corresponding form of 
the velocity potential (eqs. 3-4) and by the random walk solution illustrated in fig. 1; as noted by Weinberg & Gunn 
(1990a) the non-linear Hopf-Colc transformation ijj ^ S, in the solution of the Bernoulli eq. (1) amplifies gradients so 
that the diffusive smoothing has the greater impact where it is most needed. The lattcir point is shared by EA; however, 
the graphical illustration (fig. 1) of the EA solution for the velocity field shows that, for sufficiently large times at the 
onset of the term Vi, a second process overlaps to diffusive smoothing in determining the evolution of LSS; in fact, the 
velocity field at a point is now influenced by the value of the field at times closer to present, as shown in detail in fig. 1. 
Thus, regions with larger overdensities (corresponding to larger velocities in the surroundings) aquire a larger and larger 
role in driving the evolution of ip, and this has the effect of enhancing the disomogeneities in the density field. This is 
apparent also from the path- integral solution in eq. (20), which shows that the random diffusion term (the "kinetic" 
part of the action oc (dh/ds)'^ which corresponds to the "smoothing" mode) is now complemented by the potential term 
T^(^); of all the paths that contribute to the integral, those passing through maxima of V{^) (at evolved times) make 
the dominant contribution to the integral. Thus, the potential term in eq. (20) acts like a "selection" term (in contrast 
with the "smoothing" term driving the whole dynamics in AA) which progressively weights the more non-linear regions 
in determining the evolution of ^. It is this latter effect which drives the fragmentation of filaments into several knots (see 
figs. 2-4) observed at late times in N-body simulation and missing from the AA dynamics. Inspection of fig. 5 shows that 
such second mode in the LSS evolution begins to efficiently overlap to the diffusive smoothing already at z ^ 1. Thus, the 
non- linear gravitational evolution can be viewed as smoothing process of the initial conditions for 5 < 5, as suggested by 
previous works on AA; for larger overdensities at 2; < 1 the velocity gradient induced by small-scale overdensities overlaps 
to the smoothing mode so that particles flow along the filaments to enhance small-scale overdensities, partially breaking 
the extended structures into dense knots. A finer description of such process, includes the effect of each particle trajectory 
in the modification of the velocity field and corresponds to the term V2 in eq. (14). 

As to the investigation of cosmological problems, EA can be used to complement N-body simulations in several ways: 
first, the shorter computational time taken by EA to run allows a more extensive exploration of the parameter space in 
many astrophysical problems; second, EA can be used in the developement or testing phase of investigation techniques 
which require a large number of simulations; third, it can be used to estimate the probability for a given configuration (both 
in density and in velocity) to occur, a problem which may also require running a large number of simulations. It must be 
noticed that the above advantages are shared with other approximations like ZA and AA; however, EA allows to describe 
a wider range of overdensities (see, e.g. fig 5), thus extending the fields of investigation and allowing to address additional 
problems. A first example is constituted by the study of Lya regions and, in general, by the comparison of theoretical 
predictions with the observations concerning the intergalactic medium. Indeed, several authors have used the ZA (with 
an appropriate smoothing on initial conditions) to study the distribution of Lya column densities (Hui, Gnedin, & Zhang 
1997; Gnedin & Hui 1998); the density and velocity field derived from ZA were related to the gas density, temperature 
and composition by an independently derived equation of state. The resulting distribution of Lya column density can be 
compared with observations and a large variety of parameters (con(;erning the cosmology, the equation of state of the gas, 
the reionization epoch and the ionazing radiation) can explored through the use of the fast ZA algorithm. However, such 
approach could be applied only for density contrast S < 5 (corresponding to column densities < 10^^'^ cm~^) due to the 
break down of ZA (and also of AA as shown by fig. 5) at higher density contrast. In this context, EA could extend the 
range of such kind of investigation to larger density contrast and hence to larger column densities. In this context, such 
investigations could be further extended to include (at least partially) baryons at temperatures ^ 10"'' — 10^ K residing 
in higher density contrast 5 ^ 10, which could constitute a relevant (if not major) fraction of all existing baryons (see, 
e.g., Cen & Ostriker 1999). While a full treatement of them requires full hydrodynamical simulations including shock 
heating, preliminary studies and parameter exploration concerning the statistics of column densities of such gas could 
be performed through EA. Once a smaller set of plausible models is identified with this technique full hydrodynamical 
simulations can be run. Further examples of cosmolgical studies through EA can be constituted by the computation of 
the density distribution produced in the mildly non-linear regime cxtc^nding to density 5 10 in a variety of cosmological 
conditions; this enters many analytical or semi-analytical computations concerning the thermal and chemical state of 
the intergalactic medium (which are also being included in semi-analytic models of galaxy formations, see Benson et al. 
2001). Additional applications concern the extension to larger overdensities of a reliable velocity-density relation (widely 
investigated with the use ZA up to densities S ^4, see Nusser, Dekel, Bertshinger & Blumenthal 1991; see also Weinberg 
& Gunn 1990b), particulary used for the analysis of large-scale flows and for the inverse problem of deriving the velocity 
field corresponding to observed galaxy distribution. 

Besides complementing the N-body simulations in evolving numerically a dark matter field, the compact, anlytical form 
(eq. 14) for the evolution of the velocity field constitutes a promising way to study directly and analytically relevant 
scaling properties for the collisionless fiuid. In particular, it is known that the solution of equations similar to (14) with 
a random source term V show interesting fractal (Brax 1992) and intermittency (Gartner & Molchanov 1992) properties. 
Indeed, an approach involving a diffusion equation with a source term for the velocity field has been used by Jones (1999) 
to relate the baryonic velocity field (the one following the diffusion equation in such model) to the dark matter potential 
(the source term) which, in this approach, is given as an input. The intermittency and fractal properties of the baryonic 
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velocity field in this model (which give rise to nice scaling properties of the resulting galaxy distribution) are probably 
features which are shared by the velocity field in our model. While the investigation of such issues is more complicated in 
EA than in AA or ZA due to the more complex form of the velocity potential, it is nevertheless interesting to study the 
effects of the non-linear source term introduced by EA in the Burgers equation on the fractal and intermittency properties 
of the resulting velocity field, since this would provide useful analytical tools to characterize the growth of LSS in a more 
evolved non-linear stage than that probed by previous approximations. The investigation of such point will be the subject 
of a next paper. 

I thank S. Borgani and S. Matarrese for helpful discussions, and the referee for his constructive comments. 
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APPENDIX 

The quantity Qa = Jc^{x) ^^'^C evolves according to the following relation 

(9a+da(x) = Qa(z) + A ln^{z)da , (A-1) 

where z is the particle position at the previous time in the trajectory C{x), so that z = x — Ua{z)da. For small time 
increments da, the velocity in z is given by 

Ua(z) = u,(x) [1 - V • Ua(x)da] . (A-2) 

Substituting for z and for Ua(z) into cq. (Al), one obtains 

(3a+da(x) = ga(x) - [VQ • Ua](x) da + A ln^{x)da + o(da2) , (A-3) 

from which, after substituting u = V{—2uln^) (by definition of eq. (12) follows. 
As for eq. (13), we note that 



x = 



q-2u f V«<(x(q,a'))rfa' (A-4) 

J an 



Prom this it follows that 



A /nC(x(q, a + da)) = A Z< - 2i/ / V Z<(x(q, a')) da' - 2vV Z<(x(q, a)) da\ . (A-5) 

For small displacements — 2i/V Zn^(x(q, a)) rfa along the particle trajectory, the expansion of the argument of the r.h.s. 

yields 

A Z<(x(q, a + rfa) = A Z<(x(q, a)) - 2!/[V(A InC) ■ V Z<](x(q, a)) da . (A-6) 
whose iteration leads to eq. (13). 



